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Abstract 

We study the thermodynamical properties of a self-gravitating gas with two or more 
types of particles. Using the method of linear series of equilibria, we determine the struc- 
ture and stability of statistical equilibrium states in both microcanonical and canonical 
ensembles. We show how the critical temperature (Jeans instability) and the critical 
energy (Antonov instability) depend on the relative mass of the particles and on the di- 
mension of space. We then study the dynamical evolution of a multi-components gas of 
self-gravitating Brownian particles in the canonical ensemble. Self-similar solutions de- 
scribing the collapse below the critical temperature are obtained analytically. We find 
particle segregation, with the scaling profile of the slowest collapsing particles decaying 
with a non universal exponent that we compute perturbatively in different limits. These 
results are compared with numerical simulations of the two-species Smoluchowski-Poisson 
system. Our model of self-attracting Brownian particles also describes the chemotactic 
aggregation of a multi-species system of bacteria in biology. 

1 Introduction 

In previous papers of this series PP-|H|; we have introduced a model of self-gravitating Brownian 
particles and we studied its equilibrium and collapse properties in the framework of thermody- 
namics. In this model, the motion of the particles is described by coupled stochastic equations 
(one for each particle) involving a friction and a random force in addition to self-gravity. The 
friction and the random force mimic the influence of a thermal bath of non-gravitational origin 
imposing the temperature. The temperature of the bath measures the strength of the stochastic 
force. The self-gravitating Brownian gas model has a conceptual interest in physics because it 
represents the canonical counterpart of a Hamiltonian system of stars in Newtonian interaction. 
Therefore, it can be used to test dynamically the inequivalence of statistical ensembles which 
is generic for systems with long-range interactions. Although most astrophysical systems are 
described by the Newton equations (without dissipation), the self-gravitating Brownian gas 
model could find applications for the transport of dust particles in the solar nebula and the 
formation of planetesimals by gravitational instability 0. In this situation, the particles expe- 
rience a drag force due to the friction with the gas and a stochastic force due to turbulence. 
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Furthermore, self-gravity must be taken into account when the particles have grown sufficiently 
by sticking processes and start to feel their mutual attraction. This would be just a first step 
because other ingredients are required to improve the description of planetesimal formation. It 
has also been shown in jTUj that the process of violent relaxation for collisionless stellar systems 
exhibits similarities with the dynamics of a self-gravitating Brownian gas. In particular, the 
coarse-grained distribution function /(r, v,t) satisfies a generalized Fokker-Planck equation, 
involving an effective diffusion and an effective friction taking into account the peculiarities of 
the collisionless evolution. 

Our model of self-gravitating Brownian particles has also interest for systems that are not 
necessarily related to astrophysics. For example, in the physics of ultra cold gases, it has been 
shown recently that, using a clever configuration of lasers beams, it is possible to create an 
attractive 1/r interaction between atoms ^Tj. This leads to the fascinating possibility of re- 
producing gravitational instabilities in the laboratory. In particular, it is argued in [15| that it 
should be possible to observe the "isothermal collapse" ^3 E] of a Fermi gas cloud in thermal 
equilibrium with a bosonic "reservoir". Since the system is essentially dissipative a canonical 
description (fixed T) is required and a plausible dynamical description of the system would 
be formed by the Fokker-Planck equation coupled to the gravitational Poisson equation. In 
that case, the quantum nature of the particles (fermions) is important, and generalized Fokker- 
Planck equations, including the Pauli exclusion principle, must be considered as in (6j. On the 
other hand, as discussed in our previous papers, the collapse of the self-gravitating Brownian 
gas is analogous to the chemotactic aggregation of bacterial populations in biology. In partic- 
ular, the Smoluchowski-Poisson system which describes self-gravitating Brownian particles in 
a strong friction limit is isomorphic to a simplified version of the Keller-Segel model |J^J in 
biology, obtained in the limit of large diffusivity of the chemical JH]- The Keller-Segel model is 
a standard model in mathematical biology Due to this analogy, the results of p]-[H] have 
direct application for the chemotactic aggregation of bacterial populations. 

For all these reasons, and also in its own right, the study of the self-gravitating Brownian gas 
model P] is clearly of interest in physics. So far, most works have focused on the case of a single 
species of particles. In this paper, we extend these approaches to the case of a multi-components 
system, with particular attention devoted to the two-species model. In Sec. |21 we present 
the basic equations describing a multi-components self-gravitating Hamiltonian and Brownian 
system and show the analogies of the latter with a multi-components chemotactic system. 
We use a mean-field approach which is exact in a suitable thermodynamic limit iV„ — >■ +00, 
keeping rja = PGrnl^Na/ W^^"^ constant for each species a (see Appendix [X)). In Sec. we 
discuss the statistical equilibrium states of a two-components self-gravitating system in both 
microcanonical and canonical ensembles. Therefore, our static study applies both to ordinary 
stellar systems (galaxies, globular clusters,...) described in the microcanonical ensemble and 
Brownian systems (or bacteria) described in the canonical ensemble. We obtain the equilibrium 
density proffies and analyze their thermodynamical stability by drawing the linear series of 
equilibria (caloric curves) and using the turning point argument PEI . We show how the critical 
temperature (Jeans instability JH]) and the critical energy (Antonov instability P^ I2Uj) depend 
on the parameters fi = mi/ 1712 and x = M1/M2 where rria is the individual mass of the particles 
and Ma = Nanta the total mass of species a. In the microcanonical ensemble, we find that the 
gravothermal catastrophe is advanced, i.e. it occurs sooner with respect to the single-species 
case. In the canonical ensemble, the isothermal collapse is advanced if we add heavy particles 
in the system and delayed if we add lighter particles (keeping the total mass fixed). Exact 
analytical expressions of the critical temperature of collapse are given in dimension d = 2. An 
approximate expression is obtained for c? > 2 by using the Jeans swindle (see Appendix^. Our 
static study (Sec. EI) completes previous investigations by Taff et al. ^T] and De Vega & Siebert 
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|22j in c? = 3 and Yawn & Miller |23] in d = 1. In Sec. 01 we consider for the first time the 
dynamics of the two-species self-gravitating Brownian gas in a strong friction limit described by 
the Smoluchowski-Poisson system. We study the collapse below Tc by looking for self-similar 
solutions. Depending on the values of /x and ( = ^1/^2, where are the friction coefficients, 
we show that the collapse of one species of particles dominates the other. The invariant profile 
of the dominant species scales as p ~ as for the one-component gas PP. The collapse of the 
other particles is slaved to the collapse of the dominant species. This decouples the equations 
of motion and reduces the problem to the study of a single new dynamical equation. We show 
that this equation possesses self-similar solutions and that the scaling profile scales as p ~ r~° 
where a is a non-trivial exponent depending on p, ( and d, which leads to particle segregation. 
We determine this scaling exponent perturbatively in a large dimension limit d —>■ +00 on the 
one hand and for a weak asymmetry p — 1 and — 1 on the other hand. We also consider 
the limits p — 0, +00 or ^ — >• 0, +00. These perturbative analytical results are compared with 
the exact results obtained numerically. 



2 Analogy between self-gravitating Brownian particles 
and bacterial populations 

2.1 Self-gravitating Hamiltonian systems with different species of 
particles 

Let us consider a Hamiltonian system of X species of particles with mass mi, mx in a space 
of dimension d. Throughout the paper, the particles of species 1 are labeled from i = 1 to 
^^lA^, the particles of species 2 from {giN + 1) to g2N, and so on up to species X. The Latin 
letters i will index the N particles and the Greek letters a will index the X species. The 
particles interact via a long-range potential U{ri, r^) = J2i<j 'mirrijuijCi — rj). In this paper, 
u{ri — Tj) = —G/[{d — 2)\vi — rj\^'^~'^^] denotes the gravitational potential of interaction in d 
dimensions. This Hamiltonian system is completely defined by the equations of motion 

dvi 

dt 

-7- = V^t/(ri,...,rjv). 

dt rrii 

In kinetic theory, the coUisionless evolution of this system is governed by the Vlasov-Poisson 
system, which is valid for sufficiently "short" times. In fact, this regime can be extremely long 
in practice since the relaxation time (Chandrasekhar's time) increases almost linearly with the 
number of particles. The coUisional regime is usually described by the Landau- Poisson system 
which governs the evolution of the distribution function /(r, v, t) toward statistical equilibrium. 
For a multi-species system in ci = 3, the Landau equation reads 

^ + ^.^ + F.^ = —\^ fK^-^(m f'^-m f ^] dV 

dt ^ dr ^ ^^r dv>^ ^ l ^^^^t;- "^"^t;'- J ' 

7=1 " ^ ' 

K^"" = 27iG^- In A 6"' - ^ , (2) 
u \ ) 

where u = v — v' is the relative velocity of the particles involved in an encounter, InA = 
Jq^°° dkjk is the Coulomb factor (which must be appropriately regularized) and F = — V$ is 
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the gravitational force per unit of mass. We have also set = /^(r, v',t) assuming that the 
encounters can be treated as local (see |21I for a critical discussion of this approximation). The 
gravitational potential $(r, t) is determined by the Poisson equation 

A$ = SdGp, (3) 

with the total density p = X]a=iPa; where pa{r,t) = J /^(r, v, t)(iv is the spatial density of 
species a and /^(r, v, t) is their distribution function (/^(r, v, t) c/r (iv gives the total mass of 
particles of species a with position in (r ; r + dr) and velocity in (v ; v + rfv) at time t). The 
total distribution function is / = /«. 

The Landau-Poisson system conserves the total mass 



Ma= Padr = NaUla, (4) 



of each species and the total energy 

E = ^j fv'^drdv + ^j p<^dr = K + W, (5) 

where K is the kinetic energy and W is the potential energy. Furthermore, the Landau-Poisson 
system satisfies a H-theorem S > for the multi-components Boltzmann entropy 

X 

(6) 



S = —ks'S^ /"^^Inf^^l drdw. 



At equilibrium, S = implying that the current in the R.H.S of Eq. Q must vanish. The 
advective term in the L.H.S of this equation must also vanish, independently. These two 
conditions imply that the only stationary solution of the Landau equation Q is the Maxwell- 
Boltzmann distribution 

/„(r.v)=A<,(!^)"%-^-l*-»H. (7) 

where the inverse temperature f3 = l/ksT appears as an integration constant. Note that the 
advective term (Vlasov) is canceled out by any distribution function = /^(e) depending 
on the particle energy e = y + ^(i^) alone. The cancellation of the collision term singles out 
the Boltzmann distribution among this infinite class of distributions. The Maxwell-Boltzmann 
distribution Eq. ((Tj) represents the statistical equilibrium state of the system in a mean-field 
approximation. It can be obtained alternatively by maximizing the entropy ® at fixed energy 
and particle number (for each species). The condition of thermodynamical stability in the 
microcanonical ensemble {maximum of S at fixed E, Na) is equivalent to the linear dynamical 
stability with respect to the Landau-Poisson system |2S]- 

According to the theorem of equipartition of energy (which remains valid here), the r.m.s. 
velocity of species a decreases with mass such that 
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{v^)a = -2 = d^. (8) 

Jg-/3mc%-^^ nia 

Therefore, heavy particles have less velocity dispersion to resist gravitational attraction and 
will preferentially orbit in the inner region of the system. This leads to mass segregation, but 
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of a very different nature from the dynamical segregation that we study in Section 4. Defining 
the pressure by p = ^ / fv^dv, we get from Eq. (jH)) the local equation of state 



X 

rrir 



p = J2—kBT. (9) 



The local mass density pa of each species is obtained directly from the integration of Eq. ((Tj) 
over the velocities yielding 

p,(r) = A.e-'^"^'^*^. (10) 

The gravitational field $(r) is obtained self-consistently by substituting Eq. (fTUI) in the Poisson 
equation (jSj) and solving the resulting differential equation. 

2.2 Self-gravitating Brownian particles with different species of par- 
ticles 

The Hamiltonian system of stars presented in Sec. 12.11 is associated to the microcanonical 
ensemble (fixed energy) in statistical mechanics. We shall now introduce a model of particles in 
Newtonian interaction associated with the canonical ensemble (fixed temperature). Specifically, 
we consider a system of self-gravitating Brownian particles belonging to X different species. 
This is the generalization of the model introduced in p. This system is characterized by 
coupled stochastic equations 

-77 = Vi, 

d! 1 ^^^^ 

-T^ = V,t/(ri,...,r^) + v^Ri(t), 

at iTLi 

where is the friction coefficient, Di is the diffusion coefficient and Rj(t) the stochastic force. 
In this paper Ri(t) is a white-noise satisfying the conditions (Ri(t)) = and {Ra,i(t)Rbj(t')) = 
SabSijS(t — t') where j refers to the particles and a, h to the space coordinates. The diffusion 
coefficient and the friction coefficient are related to each other by the Einstein relation (see 
Appendix \^ 

= (12) 
rria 

where T is the thermodynamical temperature. Therefore, the temperature measures the 
strength of the stochastic force. 

In the mean-field approximation, the evolution of the system is governed by the multi- 
components Kramers equation (see Appendix El for details) 

+ V • -t:— + F ■ — — = ■ Da^::— + ^afa^ , (13) 



dt dr dv dv \ " dv 

which must be coupled consistently with the Poisson equation, using F = — V$. The Kramers- 
Poisson system conserves the total mass of each species. Since the system is dissipative, the 
energy © is not conserved and the entropy © does not increase monotonically. However, 
introducing the free energy 

F[{fa}] = E[{U}]-TS[{U}l (14) 

the Kramers- Poisson system satisfies a sort of canonical H-theorem F < 0. At equilibrium, 
F = implying that the diffusion current in Eq. (fT!?|) must vanish. The advective term must 
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also vanish. These two conditions lead to the Maxwell-Boltzmann distribution ((7j) where 1/(3 
is the temperature of the bath. This distribution represents the statistical equilibrium state of 
the system in a mean-field approximation. It can be obtained alternatively by minimizing the 
free energy ()14|) at fixed particle number (for each species). The condition of thermodynamical 
stability in the canonical ensemble {minimum of F at fixed A^^^) is equivalent to the linear 
dynamical stability with respect to the Kramers- Poisson system |25j . 

In order to simplify the problem further, we consider the strong friction limit and let 
-|-oo for each particle i. This amounts to neglecting the inertia of the particles. Instead of 
Eq. (fTTj) . we obtain a simpler system of coupled stochastic equations 



= -fi,ViU{ri, r^) + v^Ri(t), 



(15) 



where fj,i = l/rrii^i is the mobility and D[ = Di/^f = kBT/rrii^i is the diffusion coefficient in 
physical space. The mean-field Fokker-Planck equation obtained in this limit of strong friction 
is the Smoluchowski equation, which can be written for each species 



dpa 
dt 



knT. 



rrir 



(16) 



It has to be solved in conjunction with the Poisson equation Q. The passage from the Kramers 
to the Smoluchowski equation can be made rigorous by using a Chapman-Enskog expansion 
(see j2ni for details and generalizations). In the +oo limit, the distribution function can 

be written 



13 
2ti 



d/2 



Pa{r,t)e 



+ 



(17) 



where pQ,(r, t) evolves according to Eq. (fTBj) . Using Eqs. (fT^ and (fT7|) . it is possible to express 
the free energy as a functional of the spatial density of each species in the form 



p$ dr 



In 



Pa 



dr. 



up to an irrelevant additive constant. The Smoluchowski- Poisson system conserves the total 
mass of each species and decreases the free energy F < 0. At equilibrium, the density is given 
by Eq. ()1U|). The linearly dynamically stable steady states minimize the free energy F[{pq}] at 
fixed mass (for each species) pS]. 

The Kramers and Smoluchowski equations can be written 



dfa 

dt 



+ V 



dr 



+ F 



dfa 

dv 



d_ 
d\ 



(19) 



dpa 

dt 



6F 

5pa 



(20) 



where the free energy is respectively given by Eqs. (jl4j) and (jl8p . They can also be obtained from 
the linear thermodynamics of Onsager or by maximizing the rate of free energy dissipation under 
appropriate constraints j2Zl, which is the variational version of the linear thermodynamics. 



6 



2.3 Mult i- components chemotactic systems 

In previous papers, see e.g. [B], we have shown that the equations describing the dynamics 
of self-gravitating Brownian particles in a strong friction limit were isomorphic to a simplified 
version of the Keller-Segel model jTHI describing the chemotactic aggregation of bacterial pop- 
ulations. We shall propose here a simple generalization of this model to a multi-components 
system of bacteria and show the relation with the multi-components Brownian model intro- 
duced previously. Note that a more general multi-components chemotactic model has been 
proposed recently by Wolansky [2H] • We consider a system of X populations of bacteria with 
density pa, each species secreting a substance (chemical) with density Cq. The bacteria dif- 
fuse with a diffusion coefficient and they move along the (total) concentration of chemical 
c = J2a as a result of a chemotactic attraction. The chemicals, produced by the bacteria 
with a rate a, are degraded with a rate h. They also diffuse with a diffusion coefficient D' . The 
evolution of the system is described by the coupled differential equations 

^ = D«Ap,-x.V(p,Vc), (21) 
dc 

= D'Ac^ + ap^-bc^. (22) 

Like in the one-species problem ^B] , we shall consider a regime of large diffusion of the chemicals 
so that we ignore the temporal derivative in the second equation. We shall also take 6 = 0, 
assuming that there is no degradation of the chemicals. This reduces the problem to the coupled 
system 

^ = D^Ap^-XaV{paVc), (23) 
Ac = -\p. (24) 

These equations are isomorphic to the multi-components Smoluchowski-Poisson system p6|) -(Pjl 
provided that we make the identification Da = ksT / ^a^ria, Xa = ^/^a, c = — $ and A = SdG. 
Due to this analogy, the following results can be applied to the chemotactic problem in biology 
by a proper reinterpretation of the parameters. 

3 Statistical equilibrium states of a multi-components 
system of self-gravitating particles 

3.1 The thermodynamical potentials 

At a fundamental level, the Boltzmann entropy is defined by 5 = /cbIuVT where W is the 
number of microstates (complexions) associated with a given macrostate. This number W 
can be obtained by combinatorial analysis. In the continuum limit, a macro-state is specified 
by the smooth distribution function /(r,v) and the Boltzmann entropy takes the form of 
Eq. ®. Therefore, if we assume that all microstates are equiprobable for an isolated system 
at equilibrium (microcanonical ensemble), the optimal distribution function maximizes the 
Boltzmann entropy at fixed total energy and mass (for each species). Introducing Lagrange 
multipliers and writing the variational principle in the form 

X 

6S - I35E - XJMa = 0, (25) 

a=l 
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we obtain the Maxwell-Boltzmann distribution ((7j). It is important to recall at that stage that 
the Boltzmann entropy has no global maximum for self-gravitating systems. Hence, we have 
to confine the system within a restricted region of space and look for local entropy maxima. 
These metastable states are physically relevant because their lifetime increases exponentially 
with the number of particles |29j . 

On the other hand, if the system is in contact with a heat bath fixing the temperature 
(canonical ensemble), the statistical equilibrium state minimizes the free energy F = E — TS 
at fixed mass (for each species). Introducing Lagrange multipliers and writing the variational 
principle in the form 

X 

K5Mo, = 0, (26) 

a=l 

we obtain the Maxwell-Boltzmann distribution ((7j) as in the microcanonical ensemble. What we 
have done essentially is a Legendre transformation to pass from the entropy to the free energy, 
as the temperature is fixed instead of the energy. Here again, the system must be confined 
within a box and only local minima of free energy exist. 

The statistical equilibrium distribution of particles is given by Eq. (|T0|) where the gravita- 
tional potential satisfies the multi-species Boltzmann-Poisson equation 

X 

A$ = SdG J2 ^ae"^™"*. (27) 

a=l 

In the microcanonical problem (Hamiltonian systems), the inverse temperature must be related 
to the energy while in the canonical problem (Brownian systems) it is imposed by the bath (and 
the corresponding mean- field energy is interpreted as the averaged energy). Then, we can plot 
the series of equilibria (3{E). The stability of the system can be settled by the turning point 
argument jTH] as in the single-species case. Although the critical points of constrained entropy 
and constrained free energy yield the same density profiles, the stability limits (related to the 
sign of the second order variations) will differ in microcanonical and canonical ensembles. As 
these results on the inequivalence of statistical ensembles have been extensively discussed in 
the single-species case [201 we shall not go into much details here and rather focus on the 
new aspects brought by the consideration of a distribution of mass among the particles. We 
also recall that for systems with long-range interactions, the mean-field description is exact (see 
Appendix ^ so that our thermodynamical approach is rigorous. 



3.2 The two-species Emden equation 

From now on, we restrict ourselves to a system with only two species of particles with mass 
mi and m2. We assume that mi > m2 and set = mi/m2 > 1. In order to determine the 
structure of isothermal spheres, we introduce the function = m2/3($ — $o) where $o is the 
gravitational potential at r = 0. The density profile of each species can then be written as 

p,=p,{0)e-^^ , p2 = P2(0)e-^ (28) 

where pi(0) and P2(0) denote the central density. Restricting ourselves to spherically sym- 
metric solutions and introducing the notation ^ = {SdG(3m2P2{^)Y^'^f' , the Boltzmann-Poisson 
equation ()27p takes the dimensionless form 

1 d /^.-l^^=e-^ + A;.e-^ (29) 
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where A = ?T-i(0)/n2(0) is the ratio of the central numerical density rta = pa/^a of the two 
species. Equation represents the two-species Emden equation in d dimensions. It must be 
supplemented by the boundary conditions 

^(0) = ^'(0) = 0. (30) 

The one-component case is recovered for A = 0. 

The two-species Emden equation ()29|) in dimension d = 3 has been studied by Taff et al. 
|21j who plotted the density profiles and the caloric curves for different values of p. In their 
work, the ratio A of central densities is maintained fixed along the series of equilibria. We shall 
extend their study in a space of dimension d (with particular emphasis on the critical dimension 
d = 2) and consider the more physical (and more complicated) case where the ratio x = M1/M2 
of the total mass of each species (which are the conserved quantities) is kept fixed instead of 
A. This makes possible to use the caloric curve f3{E) to settle the thermodynamical stability 
of the system using the turning point argument (this is not possible when x varies along the 
series of equilibria). Furthermore, we shall obtain analytical expressions of the critical points 
(energy and temperature) as a function of fi and x- 

We shall first derive general properties of the differential equation ()29|). For ^ — > 0, an 
expansion of ipl^) in Taylor series yields 

"^^^^ " 8d{d + 2) ^ ^'^^^ 

^ ^ -[d{l + XfiY + {d + 2){l + A/i)(l + Xfi^)]e + 0{e). 



48rf2(c; + 2)(d + 4)' 

To investigate the asymptotic behavior of for ^ +00, we first perform the transformation 
t = hi^ and z = —ip + 2\ia^. In terms of z and t, the two-species Emden equation (jSHl) becomes 

^ + (rf - 2)^ = -A/ie^^e-'('^-i)* - + 2{d - 2). (32) 
dt'^ dt 

For ^ — >■ -|-oo, i.e. t — >■ -|-oo, the concentration of heavy particles, proportional to e"'^'^, goes 
to zero faster than the concentration of light particles, proportional to , so the first term in 
the R.H.S. can be neglected in a first approximation. Then, Eq. reduces to the equation 
obtained for a single type of particles. For c? > 2, it describes the damped motion of a fictitious 
particle in a potential V{z) = — 2{d — 2)z where z plays the role of position and t the role 
of time. For t — >• -|-oo, the system has reached its equilibrium position at Zq = \ia[2{d — 2)]. 
Returning to initial variables, we find that ~ 2{d — 2)/^'^ for ^ +00. Since the two-species 
Emden equation does not satisfy a homology theorem, this solution is only valid asymptotically. 
It does not form a singular solution of Eq. when A 7^ 0, contrary to the one-component 
case j2]. We note that, for d > 2, the total mass M2 ~ J^°° p2r'^~^dr of the lightest particles is 
infinite (as in the single species case) since p2 ~ r"^. However, since pi ~ r~^'^, the total mass 
of the heaviest particles is finite if 

d 

P > P3/2 = 2" (33) 

The next order correction to the asymptotic behavior of can be obtained by setting 
z = Zq + z' with z' <^ 1 and keeping only terms that are linear in z'. This yields 

'^^ + (d-2)^ + 2(d- 2)z' = -A/i2^(rf - 2)^e-2(^-i)*. (34) 
dt'' dt 
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Figure 1: The dimensionless density profiles pi{C,) = Xjie^^"^ and P2(0 = for = 5 and for 
A = 1 in = 1 (Fig. QJa) and c? = 3 (Fig. ^b). The dashed hne represents the density of the 
one component system. 

This differential equation can be solved analytically. The discriminant associated to the homo- 
geneous equation exhibits two critical dimensions d = 2 and d = 10 j2]. For 2 < c? < 10, we 
have for +oo, 



2{d-2) 



A 



■ cos 



v/(c/-2)(10-rf) 



(35) 



\^i2^'-\d - 2^^ 



Mf-2(M-1) 



2(/i-l)2-(d-2)(/i-2) 



where A and S are integration constants. The density profile ()35p intersects the asymp- 
totic solution 2{d — 2)/^^ at points that asymptotically increase geometrically in the ratio 

^ . ^2-K/^{d-2){W-d) _ 



H > /i5/4 



d + 2 



(36) 



the last term in Eq. ()35p can be neglected for sufficiently large ^ and there is an infinite number 
of intersections. For /i < /i5/4, there is only a finite number of intersections. For d > 10, we 
have for ^ +oo. 



2{d-2) 



V(d-2)(d-10) 



B 



_ V(d-2)(d-10) 



(37) 



A/i2'^-i(rf- 2)^^-2(m-i) 
"2(/x-l)2-(d-2)(/i-2)J 

There is no intersection with the asymptotic solution. For d < 2, the density profile of the 
lightest particles decreases as ~ Ce~^'^^ and for c? = 2 as e""^ ~ AC,^^. The normalized 
density profiles are plotted in Fig. [T] in c? = 1 and d = 3. The case = 2 is postponed to 
Sec. Eli 



3.3 The Milne variables 

Since the multi-species Emden equation does not satisfy a homology theorem, it cannot be 
transformed into a first order differential equation as in the one-species case. However, the use 
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of the Milne variables is still useful to analyze the phase portrait of the equation. In the general 
case, they are defined by 

_ d In M (r) _ Xfie-f"^ + e-^ _ d\Yip{r) _ ( Xpe-^"^ + e-^\ 

dlnr iJ' ' ''""^h^"^'^ ^ Ae-z^^ + e-^ J' ^^^^ 

where we used the integrated density M(r) = pSdr'^~^ dr and the total pressure p{r) = 
{pi/rrii + p2/m2)kBT . Taking the logarithmic derivatives of u and v with respect to ^ and 
introducing the notation w = {Xe~^^^ + e~'^){Xp^e~^'^ + e~^)/{Xpe~^^^ + e~'^)^, we get 

1 du 1 , 
-- = -id-u-vw), (39) 

= i[« + .;(l-^)-(rf-2)]. (40) 

The single species case is recovered for A = and w = 1. Taking the ratio of the above 
equations, we obtain 

u dv u + v(l ~ w) — (d — 2) , , 

-— = ^- ^ '-. 41 

V du d — u — vw 

The solution curve in the [u, v) plane is plotted in Fig. |21 The curve is parameterized by C,- 
It starts from the point {u, v) = {d, 0) with a slope 

'dv\ ^_d + 2 (1 + A/i)^ 

duj^ d (l + A)[l + Ap(/i(l + A/i) + l)]' ^ ^ 

corresponding to ^ = 0. For d > 2 and ^ +oo, the curve converges to the limit point {d — 2, 2) 
which corresponds to the asymptotic behavior e"'^ ~ 2{d — 2)/^'^. Contrary to the single-species 
case, the (m, f) curve can make loops before spiraling around the limit point. These loops are a 
signature of a multi-components system : using Eq. (PT|) . the points of horizontal and vertical 
tangent are defined respectively by m -|- v{l — w) = d — 2 and u + vw = d. Due to the term w, 
new solutions of these equations arise with respect to the single-species case and they create 
loops. For d > 10, the [u, v) curve reaches the limit point without spiraling but still makes 
loops for the reason previously mentioned. For d < 2, the curve tends monotonically to (0, +oo) 
for ^ — > +00 as in the single species case. The two-dimensional case is discussed in Sec. 13.51 



3.4 The thermodynamical parameters 

As indicated previously, isothermal self-gravitating systems have infinite mass. We shall over- 
come this problem by confining the system within a spherical box of radius R (Antonov model). 
Physically, the box delimits the region of space where thermodynamical arguments can be ap- 
plied. In the biological problem (chemotaxis), the box represents the natural boundary of the 
domain in which the bacteria live. For bounded isothermal systems, the solution of Eq. (|29p is 
terminated by the box at a normalized radius given by a = {SdGm2pp2{0)Y^'^R. We shall now 
determine the temperature and the energy corresponding to the configuration indexed by a. 
Using the Poisson equation (jSj), we write the Gauss theorem 

GM = G I pdv = SdG r pr"-^ dr = [""^f r''''^] dr = ( r''-'^] . (43) 



^0 



dr \ dr J \ dr J 



Introducing the dimensionless variables defined previously, we find that the normalized inverse 
temperature is given by 

pGMm2 . ,,,, 

^ = \d-2 = (")• (44) 
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Figure 2: The solution of the two-species Emden equation in the {u,v) plane for (i = 1 (left, 
Fig. |21a) and in d = 3 (right, Fig. |21b). The single-species case is represented by the dashed 
line. 



The calculation of the energy E = K + W is a little more intricate. The kinetic energy is 
given by 

K^liN,^ N.) k,T ^UMl^Ml) ^iMlU^ i) u,T. (45) 

2 2 \mi J 2 1712 / 

Using M = M2{x + 1) and Eq. ()44|) . the normalized kinetic energy can be written 

GM2 2a^'(a)/i(x + !)■ 
For d ^ 2, the expression of the potential energy can be deduced from the Virial theorem 

2K +{d-2)W = dVdR'^p{R), (47) 

where Vd = Sd/d is the volume of a d-dimensional sphere with unit radius (and Sd is the 
surface of a rf- dimensional unit sphere) [3]. Using p{R) = {pi{R)/mi + p2{R)/m2)kBT and the 
expressions fOH|l of the density, we directly obtain 

= I (Ae-'^^^") + e-^(")) + - (48) 

GM^ {d-2)^'^{a)^ ^ {d-2)a^'{a)fi{x + l)' 

Adding Eqs. pUj) and PHj) . we find that the total normalized energy is 

A = = ^ ( ^ + /^ ^ (Ae-^'^^") + e-^(")) (49) 

GM^ 2a^'{a)\d-2j p{x+^) {d - 2)^'\a) ^ >' ^ ' 

Note that an alternative expression of the potential energy, valid also for d = 2, can be obtained 
along the following lines. Starting from the expression 

W =1 / p$c/r, (50) 



2 _ 

and introducing the dimensionless variables defined previously, we get 
WR'^-^ 1 



GM2 2a'^V'^(tt) Jo 



{Xpe-^^ + e-^){il: + il:o)e-'d^, (51) 
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Figure 3: Series of equilibria (caloric curves) for a two-components isothermal gas in = 1 
(left, Fig. |21a) and in = 3 (right, Fig. |21b). We plot the inverse normalized temperature 
f] = GMm2l3 / R^'"^ as a function of the normalized energy A = —ER'^~'^/GM'^ for several 
values of the total mass ratio x = M1/M2. The dashed curve represents the one component 
case, i.e. x = 0. 



where ipo = w7.2/5$(0) represents the normalized central gravitational potential. It is determined 
by the relation ip{a) = ma/? ($(i?) - $(0)) with <l>(i?) = -GM/[{d - 2)R'^"^] for d ^ 2. This 
yields 

v,„ = -(i;;^^+^H). (52) 

Equation (j^Tj) remains valid for d = 2 but in that case, $(-R) = so that tpQ = —ipla). The 
corresponding expression of the total normalized energy is now 

Equations and define a series of equilibria P{E) parameterized by the value of the 
normalized radius a, or equivalently by the density contrast TZ = p2{0)/p2{R) = e'^^°'\ Along 
this series of equilibria, we can either fix the ratio of central densities A or the ratio of total 
mass X- These two parameters are related to each other by 

Ml ^ XpJ.^e-'e-^'^d^ .... 
^~ M2 j;^ ^''-'e-^ d^ ■ ^ ' 

In the framework of statistical mechanics, it is more relevant to fix x along the series of equilibria 
since the total mass of each species is a conserved quantity. Furthermore, it is only under this 
condition that the turning point argument can be used to settle the stability of the system. 
Therefore, in the foregoing equations, A must be viewed as an implicit function of a given by 

. (55) 

Then, for given a, the two-species Emden equation must be solved by an iterative procedure 
in order to satisfy the constraint (|^. 

Figure El displays an ensemble of caloric curves in d = 1 and in c? = 3 for different values of x 
(at fixed /i). In = 1, the curves are monotonic and the system is always stable. In d = 3, the 
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Figure 4: Evolution of the critical normalized inverse temperature rjc (Jeans temperature) and 
the critical normalized energy (Antonov energy) as a function of x for d = 3 and /i = 5. We 
plot with a dashed-dotted line the critical temperature rjj obtained by using the Jeans swindle 
(see Appendix |B] for details). Note this "naive" prediction provides a reasonable fit of the exact 
critical temperature rjc- 

curves present turning points at which a mode of stability is lost depending on the ensemble 
considered (a vertical tangent corresponds to a loss of micro canonical stability and a horizontal 
tangent to a loss of canonical stability). These results have been discussed in detail for the one- 
species case, see e.g. [13], and will not be repeated here. We shall just discuss how the critical 
points (beyond which no equilibrium state exists) depend on the relative mass of the particles. 
First, consider the canonical ensemble in which the control parameter is the normalized inverse 
temperature r]. For r] > rjc^fi, x), the system undergoes an "isothermal collapse". For x = we 
obviously recover the value rjc — 2.52 of the single species case. As the mass ratio x increases 
(at fixed total mass M and fi > 1), rjc decreases (T^ increases) up to rjc = 2.52//i obtained 
for X ~^ +00. In the microcanonical ensemble, the control parameter is the normalized energy 
A. For A > Ac(/U, x), the system undergoes a "gravothermal catastrophe". For x = and 
X — ^ +00, we recover the single-species value Ac ~ 0.335. Between these two extreme values, 
Ac passes by a minimum (Ac)mm(/^)- These results are illustrated in Fig. IHwhere rjc and Ac are 
plotted as a function of x for a given value of /i. The value of the minimum of the normalized 
energy (Ac)mm(/^) seems to behave linearly with fi (except for /i — 1) as illustrated in Fig. El 
If we take the particles of mass m2 as a reference, we conclude that the onset of isothermal 
collapse in the canonical ensemble is advanced when heavier particles mi > 1712 are added to 
the system (keeping the total mass M fixed). It is delayed if lighter masses (yU < 1) are added. 
On the other hand, the onset of the gravothermal catastrophe in the microcanonical ensemble 
is always advanced in a multi-species system (with respect to the single species case), whatever 
the mass of the particles added (keeping the total mass M fixed). 

Some analytical results can be obtained for a +00 and d > 2. This corresponds to the 
configurations located near the limit point in Fig. Elb, at the center of the spiral. In that case, 
it will be shown a posteriori that A(a) diverges for a fixed x- Accordingly, we can neglect the 
term e~'^ in the Emden equation (|29|) which reduces to 




(56) 
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Figure 5: Evolution of the minimum Antonov energy (Ac)mm as a function of fi for d = 3. In the 
range considered, it decreases approximately linearly as — 0.19719yU. Note that the minimum 
Antonov energy becomes positive for ^ ~ 2.627 and x — 0.68. Furthermore, the value of x for 
which Ac is minimum is always close to 0.7 (except for /i — > 1). This is probably related to the 
fact that (Ac)m,m(/^) is almost linear in this range. 



This approximation is valid for ^ < where ^' is such that A/xe ~ e ''I'^^'K If we introduce 

a new potential depending on C = VXjj,^ through the defining relation 

^(0 = (c) = (yX/ie) , (57) 

then Eq. ()56|) takes the form of the ordinary Emden equation 

1 d f^^^^de\^_^ ^^g^ 



Using the behavior B ~ 2ln( — ln[2(c? — 2)] for large (, we obtain the following behavior of 
ip{^) in the range 1 <^ ^ < 

e-* ^ PM!. (59) 

We shall find a posteriori that ^' ~ a so that the range of validity of this behavior is huge in 
the limit a +oo. This scaling in ^^'^^'^ contrasts from the scaling in obtained in the 
limit +00 for fixed A. The validity of Eq. ()59|) is confirmed in Fig. IHl where we plot the 
normalized density profile of a bounded isothermal system for a large value of a. 

Using Eqs. and (j^ . we can investigate the asymptotic behavior of X{a) for a +oo. 
We have to estimate the two integrals 



h{a) = / ^'-'e-f^^dC, (60) 
h{a) = re-'e-Ui, (61) 
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Figure 6: Dimensionless density profile P2{C,) = of tlie liglitest particles enclosed within a 
box for (i = 3, yU = 5, X = 1/9 aiid a = 5000. In the limit a +oo, the profile decays as ^"^/^^ 
for 1 <^ ^ < ^' ~ a. This can be contrasted to the decay for ^ +oo in an open system 
with fixed A (see Fig. [T}. 



for large values of a. We check that for ci > 2 and fi > 1, the integrals (extended to +oo) do 
not converge. Therefore, both terms in the decomposition 

/2(a) = ^^-^e-^d^+ [ i'^'^e-^di (62) 



Jo J^' 

behave as a power law with the same exponent (the second integral is not negligible with respect 
to the first). We can obtain the asymptotic behavior of the first integral by using the analytical 
expression (j59|) of ip in the range 1 -C ^ < However, since we do not know the expression of 
for ^' < ^ < a, we cannot compute the second integral. Thus we can get the exponent of the 
power law divergence of Ii{a) but not the prefactor. 

Evaluating Eqs. (pH) and (jUT|) with Eq. we obtain 

. «a)~A-.ig;. (63) 

where, for the reasons explained previously, the prefactors are not known. The asymptotic 
behavior of A (a) is now obtained by substituting Eq. in Eq. This yields 

X{a) ^ Ka^'^^'-^\ (64) 

As yU > 1, the numerical density ratio A(a) always diverges for a +00 and d > 2. This 
justifies our initial assumptions. If we now insert Eq. ()64|1 in ()63|1 . we find that /2(a) diverges 
as /2(a) ~ a*^"^. On the other hand, Ii{a) behaves as Ii{a) ~ a''"^^ which diverges for 
A* < /^3/2 = d/2 and tends to zero for /i > /i3/2. In Fig. [7|we plot the ratio of central numerical 
densities A as a function of t] for different values of x- 



3.5 The two-dimensional case 

The dimension ci = 2 is a critical dimension for self-gravitating systems |2] . It is also the relevant 
dimension for the biological problem of chemotaxis, since bacterial colonies usually live on a 
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Figure 7: The A (77) curves for several values of x = M1/M2 for d = 3. These curves are 
parameterized by a. 



plate. Therefore, the dimension d = 2 requires a particular attention. The two-dimensional 
Emden equation reads 

The density profile behaves asymptotically as e ~ AC^ ^ where A and 5 are constants. For 
A = (single-species case), Eq. (jU^ can be solved analytically and we get the exponent 6 = 4. 
For other values of A, we have 5 7^ 4. Some density profiles are plotted in Fig. |H1 The phase 
portrait of the Emden equation ()29p in the Milne plane is shown in Fig. El On the other hand, 
the thermodynamical parameters t] and A are given by 

ri = GMm2p = aij'{a), (66) 



Note that the normalized temperature and the normalized energy do not depend on R. This 
is a consequence of the logarithmic form of the gravitational potential in c? = 2. An ensemble 
of caloric curves are plotted in Fig. E3 As in the previous section, the value of x is fixed 
along a series of equilibria so that A(a) is determined by an iterative procedure. In continuity 
with the single species case, the caloric curves form a plateau for A — > +00. Thus, there 
exists a critical inverse temperature //c(X; A^) above which no equilibrium state is possible in 
the canonical ensemble (by contrast, there is no critical energy in c? = 2 in the microcanonical 
ensemble). The critical temperature rjclx, f^) has two expressions depending on whether > 2 
or /i < 2 as we now show. 

We first consider the situation in which, at T = Tc, the two profiles form a Dirac peak 
at r = 0. The parameters corresponding to this situation will be found a posteriori. In this 
case, the critical temperature can be obtained from the Virial theorem as in the single-species 
problem . We start from the general relation valid in c? = 2 : 

2K-^ = 2piR)V, (68) 
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Figure 8: Dimensionless density profiles pi{C,) = \pe~^'^ and P2(0 — dimension d = 2 

for A = 1 and for /i = 2 and fi = 5. For comparison, the dashed line represents the density of 
the singe-component system with slope —4. The slope of the profile ~ depends on A 
and /i. 



> 




Figure 9: The solution of Eq. (jSHI) in the (m, v) plane, where u and v are defined by Eq. pKjl for 
the two-components system at c? = 2. The single-species case is represented by the dashed hne. 
All the curves start at {u,v) = {d,0) with the initial slope given by Eq. PT|). For ^ — > +oo, 
they tend to the terminal point (0, 6). 
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Figure 10: An ensemble of caloric curves for different values of x = M1/M2 in the two- 
dimensional case. The dashed curve represents the one-component case. 

with V = TiR^. If the densities are concentrated in a Dirac peak at r = 0, the total pressure at 
the edge of the box vanishes. Equations ()45|) and (jUHj) directly lead to the result 



GM2 ^ ^ 

ksT. = -j^. (69) 



The critical normalized temperature is 



For X = 0, we recover the critical inverse temperature rj^ = ^ obtained for the single-species 
case (2]. 

It will be shown that the above regime, called regime (I), corresponds to the case where the 
function \{q) converges for a —>■ +00. We now consider the regime (II) where A(a) diverges 
so that Eq. (jU^ reduces to Eq. with d = 2. With the change of variables of Eq. (j^Tj), we 
obtain the classical Emden equation (j58|) . In d = 2, it can be solved analytically and, returning 
to original variables, we get 

m = lin(i + ^e). (71) 

This analytical expression provides a good approximation of the solution for all values of < ^' 
where ^' is such that 

e-^^^^^d^ ^ Xfi / e'^'^^^^^d^. (72) 
'0 Jo 

Substituting Eq. (f7T|) into Eq. (f7^ and using Eq. (f7i|) . we find that ^' ~ a so the range of 
validity of the analytical expression is huge, as checked numerically. The density profile of 
the lightest particles decreases as ~ C"^'''^ cind the density profile of the heaviest particles 
decreases as e~^^ ~ These heaviest particles form a Dirac peak as a ^ +00. Furthermore, 
the density \^e~^^ of species 1 becomes smaller than the density of species 2 for ^ > ^" ~ 
q,-(m-2) /(4(/x-i)) _^ g_ order to determine the asymptotic behavior of \{q) for a +00, 
we have to estimate the integrals /i(a) and h^a) defined by Eqs. (jHIH) and (jHH). We will 
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Figure 11: The critical temperature rjc is plotted versus x (Fig- CHa) and /i (Fig. [TTJb). The 
solid lines represent the numerical results. They are in excellent agreement with the theoretical 
results (I76j) and (j7U|) in the two regimes. The transition between these regimes appears when 
X = X* = l/(/^ ~ 2) in Fig. [TTJa and when /i = /i* = 2 + l/xin Fig. ITTlb. 

consistently show that the assumptions made in regime (II) are only valid for fi > 2. Then, 
it is easy to show that the integral /i extended to +oo is convergent while the integral I2 is 
divergent. Therefore, by using the profile (fTTjl to evaluate (pljl and (pT|) . we obtain the exact 
asymptotic expression of Ji while we only obtain the correct exponent of a in I2 but not the 
prefactor. Indeed, in the calculation of Ji the second integral in the decomposition (j62j) is 
negligible while in the calculation of I2 it is of the same order as the first. Then, we obtain 

A(«)~A ' /2(a) ~ ir2A-'/'^a'(^-')/^ (73) 

and, using Eq. we get 

A(a.)~(^^)*a"-. (74) 

for a —>■ +00. The assumption that A(a) diverges is only consistent with /i > 2. Note that 
from Eq. (ffljl and the value of Ji, we obtain the exact result I2 — > 4/x/i. We also note that 
/i — > for a — > +00. 

We are now able to obtain the critical inverse temperature rjc in regime (II). In this regime, 
the heavy particles form a Dirac peak at r = for T = while the light particles extend in the 
whole box. Since their density is non-zero on the edge of the box, we cannot use the reasoning 
valid in regime (I). However, the normalized temperature can be written in the form 

T] = (3GMm2 = 2nl3Gm2 + ^) ^ pi{r)rdr = Xfi + ^) ^ e"^'^^ C^^) 
where the last integral is precisely Ii{a). Using Eq. (fTSj) for a +00, we get 

Vc = -(l + -). (76) 

We note the fortunate cancellation of A which allows one to obtain this exact result without 
detailed knowledge of K2. We now have two different expressions of the critical normalized 
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Figure 12: The normalized density profile p2 = e~'^ is plotted for /i = 5 > 2 and for two 
different values of x situated from both sides of the critical value x* = l/(/^ ^ 2) = 1/3. For 
X = 1/9 < X* (regime I), the asymptotic slope of the profile is —rjc = —3.6799 given by Eq. (fTUj). 
For X = ^ > X* (regime II), the asymptotic slope in the range 1 <^ ^ ^ ^' ~ a is — 4//i = —0.8, 
see Eq. (jTT)). The final asymptotic slope in the range ^' < ^ < a is — ?7c = —8/5 given by 
Eq. d^. 



inverse temperature rjc, resp. Eqs. flTOI) and (f7n|) . We find that the crossover between the two 
regimes is obtained for 

X* = ^17 , /i* = 2 + -. (77) 

Applying the Virial theorem in regime (II), we can determine the exact expression of the 
normalized density of species 2 on the box. Indeed, using Eqs. and (f7B|). we obtain 

^ (78) 

for a —>■ +00. This implies another necessary condition to be satisfied in regime (II), namely 
X > X*- 

In conclusion, regime (I) corresponds to /i < 2 and (/i > 2 and x < X*) ! i^i that case, \{a) 
converges and the critical temperature is given by Eq. (fTUj) . Regime (II) corresponds to (/i > 2 
and X > X*) j ill that case, A (a) diverges and the critical temperature is given by Eq. (fTUI) . 
Equivalently, for a given x, if yU < yU* (regime I) the critical temperature is given by Eq. (f70|) 
while if /i > (regime II) it is given by Eq. (ffn|) . Figure ITTl clearlv exhibits the cross-over of 
these two different regimes. The theoretical predictions (ffn|) . (ffF)|l of the critical temperature 
are perfectly consistent with the numerical results. We plot the normalized density p2 = in 
Fig. El In the regime (I) where A(a) converges, we have the asymptotic behavior ~ 
for ^ — * +00. Using the fact that rjc for +oo, we find that S = rjc where rjc is 

explicitly given by Eq. (j7(Jj) . In particular, 6 = rjc = 4: for the single species case. In the regime 
(II) where \{a) diverges, we have the asymptotic behavior ~ A^~^^^ for 1 -C ^ < ^' and 
a — i> +00. The final asymptotic slope in the range ^' < ^ < a is —rjc given by Eq. (fTBj). The 
numerical results are fully compatible with these predicted values. 
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We now address the determination of the critical temperature m d = 2 for a system with 
more than two types of particles. If all the species collapse on a Dirac peak at T = Tc as 
in regime (I) discussed previously, the critical temperature is again given by Eq. (jU^ with 

= Yla-^a and M = 'Y^^Na'ma- The critical normalized inverse temperature 77 = (3GMmx 



IS 



ric = -j^. (79) 

Alternatively, we can consider the case where Ai(a) diverges for a +00, as in regime (11) 
discussed previously. To apply the same approximations as before, we need to have Ai ^ \i 
for i = 2, ...,X where Aj = ni{0)/nx{0). Repeating the steps described previously in this more 
general situation, we find that Ai — > +00 if /ii > 2 and Ai ^ Aj if fii > 2/ij. We shall assume 
that these conditions are fulfilled (i.e. mi > 2m2). In that case Ai ~ a^f^'^~'^^ and we find by an 
approach similar to that described previously that 

= -rr^. 80 

4 Collapse of a multi-components system 

4.1 Self-similar solutions of the two- components Smoluchowski-Pois- 
son system 

We now consider the dynamics of a system of self-gravitating Brownian particles. We restrict 
ourselves to the case of only two types of mass mi and m2, as we shall see that the general 
case of a discrete spectrum of particles is a simple generalization of this problem. We also 
restrict our analysis to a spatial dimension d > 2. The dimension d = 2 is critical and deserves 
a particular treatment (see 2j for the single species case). As in our previous works, we 
consider a limit of strong friction ^ — > +00 so that the dynamical equations reduce to the two- 
species Smoluchowski-Poisson system (fTBjl . We also restrict ourselves to spherically symmetric 
solutions. By introducing dimensionless variables, we can set ks = G = R = M = mi = = 1 
without loss of generality. Then, the problem depends only on the asymmetry parameters 
/i = mi/m2 = l/m2 and ( = ^1/^2 = 1/^2 and on the temperature T = l/{rj^). With these 
conventions, the dynamical equations can be written 



dp 



dt 
dp2 
dt 



V- (TVpi + piV$), 
CV- (T/iVp2 + P2V$), 



A$ = Sdp. (82) 

We shall impose a vanishing flux across the surface of the confining sphere. Therefore, the 
boundary conditions are 



dr 2-d (g3^ 

T^(l, t) + pi(l, t) = 0, Tp^(l, t) + p2(l, t) = 0. 
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Using the Gauss theorem, we can rewrite the Smoluchowski-Poisson system (j8ip - (j8'2p in the 
form of two integrodifferential equations 



dpi 


1 d 


dt 




dp2 


C d 


dt 


/yd ]_ ^^'y 



T 



dpi 
dr 



Pi 



rj. dp2 p2 

Tp— h 



dr 



-1 



Sdp{r')r"^-^ dr' 

r 

Sdp{ry^-^ dr' 



B4) 



The Smoluchowski-Poisson system (|H^ is also equivalent to a set of two coupled differential 
equations 



dMi 
dM2 



T 



d^Mi 1 - d dMi 



dr"^ 



+ 



dr 



+ 



Ml + M2 dMi 



Tp 



92 M2 l-d dM2 \ Ml + M2 dM2 



+ 



dr 



+ 



dr 



(85) 



for the quantities 



M„(r,t) 



p^{r',t)Sy''-'dr' 



which give the mass of species a = 1,2 within a sphere of radius r. In terms of these variables, 
the boundary conditions take the form 



M„(0,t) = 



Mi(l,t) 



X 



1 + X 



M2(l,t) 



1 



1 + X 



^7) 



Note that we shall restrict ourselves to the pre-collapse regime, so that we do not consider the 
possibility that a Dirac peak forms at r = 0. A Dirac peak forms in the post-collapse regime 
for (i > 2 and in = 2 (see ^E] in the single species case). It will be more convenient to work 
in terms of the functions Sa{r,t) = Ma{r,t)/r'^ which have the dimension of a density. They 
satisfy 



dsi{r,t) 

dt 
ds2{r,t) 

dt 



T 



d'^si d + l dsi 



+ 



Tp 



d'^S2 d+l ds2 



dsi , 
+ {si + S2) ( r— — h dsi 



Q^2 



+ 



r dr 



dr 



I X , ds2 , 

\si + S2) I r— + ds2 
dr 



We look for self-similar solutions of the form 



si{r,t) = po{t)Si 



S2{r,t) = po''^{t)S2 



ro{t) 



(89) 



where po{t) represents the typical central density of species 1 and ro(t) is the typical core radius 
(of the two species) defined by 

Porl = T. (90) 

On physical grounds, we expect that the total density should scale as in the single-species 
case because, on a coarse-grained scale the fine structure of the mass distribution should not 
matter (except for a continuous spectrum of mass going from [0,-|-oo[ with peculiar behavior 
at the extremes, which is not the case here). Therefore, either the two profiles scale the same 
manner or one dominates the other. Now, by solving numerically the scaling equation coming 
from Eqs. (j88|) -(|89| ) . we have found that the problem does not admit any physical solution 
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with a = 2. Hence one species will dominate the other. We define species 1 as the one that 
dominates the dynamics. This choice imposes a < 2 for the other species. We will give later 
the conditions on /i and C, for which this requirement is satisfied. Inserting Eq. ()89|) in Eq. ()88p 
and using the notation x = r/ro(t), the equation for Si(r, t) is transformed into 

+ (poSi{x) + Pq^'^S2{x)) {xpoS[{x) +dpoSi{x)). 



For sufficiently high densities, we can neglect the sub-dominant term Pq^'^S2{x) in the above 
equation. Then, Eq. reduces to 



which coincides with the equation obtained in the single-species case |2j. Setting p^^dpo/dt = 2, 
we find that 

Po{t) = 2 (tcoii - ty' . (92) 

Thus, the central density diverges in a finite time tcou- Furthermore, the differential equation 
for the invariant profile can be solved analytically j2] and we get Si{x) = So{x) where 

Using the preceding results, the differential equation determining the invariant profile of species 
2 is given by the linear second order differential equation 



CpS'^ix) + 



c(^^ + x5o(x))-x 



S'2{x) + {dCSo{x)-a)S2{x)=0. (94) 



For X +00, we have the asymptotic behavior 

S2ix) ~ x"". (95) 

Equation (j9H) can be numerically solved for any couple (/i, C)- As this equation has been 
obtained under the assumption that the exponent a < 2, we define a critical ratio of friction 
coefficients Cc(/^) corresponding to the limit of validity of this hypothesis, i.e. a{fiXc) = 2 
(similarly, we define pdC) such that a{pc,C) = 2)- 

In the case of a discrete spectrum of particle masses and friction coefficients, the above 
calculations can be repeated. One obtains an equation identical to Eq. for each type of 
sub-dominant particles, for which the analysis that we present below has to be applied. 

The critical ratio Cc(/^) is plotted in Fig. ^| The value Cc(0) can be obtained analytically. 
Inserting p = and a = 2 in Eq. (jMj) . we obtain 

xiCSoix) - l)S',{x) + (dCSoix) - 2)S2{x) = 0. (96) 

This equation must have a solution for any value of x. Taking a; = 0, we find the necessary 
condition 

C.(0) = (97) 
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Figure 13: The critical ratio Cc as a function of /i in c? = 3 is plotted in log-log scale. This 
function starts at Cc(0) = 1/6 = {d—2)/2d and diverges at yU = 2. Below the critical line, species 
1 dominates the collapse and above the critical line, species 2 dominates. In that case, our study 
can still be used with the transformation {(,fi) (l/C; 1//^)- This gives a corresponding point 
located below the dashed curve, corresponding to the function l/(c{i/ fi). We show an example 
illustrating this transformation when ( > The (x) symbol represents the point (0.65,0.41) 
and the (+) symbols represent (0.65,0.61) and (1/0.65,1/0.61) respectively. 



Inserting this result in Eq. (|^. we find that the scahng profile 5*2 (x) is 

^'<^' = (d-2)V^ + .^ - 

where A is an integration constant. The assumption that species 1 dominates the collapse is 
valid for ( < Cc(/^)- Above this critical line, the role played by the two species is swapped, 
and species 2 dominates. We can return to the studied situation by simply changing the index 
1^2. Therefore, if (/x, () belongs to region 2 in Fig. ^1 this transformation leads to study 
the case (1//U, which belongs to a sub-part of region 1, below the dashed line. 

We have numerically studied this inversion in Fig. We first start with a value of (/i, () 
below the critical line. Specifically, we take /i = 0.65 (leading to Cc — 0.51) and A = ( — (c = 
—0.1 (point x). In that case, species 1 dominates the collapse : its profile decreases as pi ~ 
while the profile of species 2 decreases as p2 ~ with a = 1.85337 < 2 determined by solving 
numerically the scaling equation with (/x, C)- We then increase the value of ( above the 
critical line, at the same distance A = ( — (c = +0.1 (point +). In that case, we have a reversal 
of population. It is now species 2 that dominates the collapse : its profile decreases as p2 ~ 
while the profile of species 1 scales as pi ~ r~" . To get the value of a' from our study, we 
set 2^1 and 1 —>■ II. We are now in the situation where species I dominates. Due to this 
transformation, the new parameters are p = mj/mji = 1/p and ( = (i/Cii = 1/C- Then, 
Si = S2 is given by Eq. (jHSI) and Su = Si is solution of Eq. (jMj) with p and (. The numerical 
solution of this scaling equation gives a' = 1.74238. We note that a' 7^ a so that the slope of 
the function a(A) — 2 is discontinuous as A — >• 0. 

For p ^ 2, the critical ratio Cc(/^) diverges. Therefore, for p > 2, species 1 always dominates 
the collapse whatever the value of It is possible to show the signature of this phenomenon 
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Figure 14: We plot the scaling profiles Si and S2 for /i = 0.65 and for different values of (. We 
take C = 0.41 < Cc = 0.511070 (dashed lines) and ( = 0.61 > Q (solid lines). This corresponds 
to the points marked x and + in Fig.^1 For ( > (c, the exponent a' is obtained from Eq. ()94j) 
using the equivalent point (l//i, 1/C). 



analytically. Assuming a = 2 and Cc +00, Eq. reduces to 

S'^+(^ + ^)s', + ^S, = 0. (99) 

\ X jl J fi 

For large x, the profile 5*2(0;) should decay as which immediately implies that n = 2. 
Equation (j^J^ for 6*2 (x) can then be solved in terms of hypergeometric functions. On the 
other hand, considering a perturbation expansion d +00 (see Sec. 14. 2p . we can obtain the 
analytical expression 

Cc(/i) = ^, (d^+oc). (100) 

We note that, in this limit, Cc(0) = 1/2 in agreement with the exact result (jU7|) . We can also 
obtain an approximate expression of the profile 5*2 (x) for d +00 (see Sec. 14. 2p . 



4.2 Perturbation expansion for d +00 

We shall first obtain the expression of the scaling exponent a for d —>■ +00. We shall see 
that the resulting expression applied for d = 3 already provides a good approximation of 
the exact solution. We use a method similar to that developed in j2] in a slightly different 
context. Equation (j94p can be formally written as a first order differential equation (writing 
5*2 = {32/52)82) depending on x, 6*0 and 82/82, 



8^ 
82 



d(8o{x) 



a 



X 



c 



'x8,{x) + J^ + ^^ 



:ior 



The term dC8Q{x) — a vanishes for a particular x, noted Xq 7^ 0, whereas the ratio 82/ 82 cannot 
vanish. This implies that the denominator in Eq. ()101|) must be equal to zero for x = xq. Using 
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Eq. (jnSl), is explicitly given by 



'd 



a 



1 



2. 



The condition that the denominator vanishes for this value can be written 



+ XqSo{xo) 



0. 



(102) 



(103) 



In the sequel, it is more convenient to work with the variable u = x'^. In terms of this variable, 
Eqs. (Uni and dTU!^ become 



Mo = ( — - 1 1 + 2 



a 



Mo - C 



UoSo{uo) 



(104) 
(105) 



Using Eqs. (jHSI) and (fTUijl . Eq. (fTUKIl can be rewritten 
'4C 



- 1 - C/i + 2 



l-C l+/^ + yUC? 



Va /Si 



(mo) a 
(mo) 2 



0. 



^5^(M0) ' d> - 

In the limit d +oo, keeping only the dominant terms in the above equation, we obtain 
+ 1 — 4(/a = from which we derive the zeroth order expression of a 

a = (107) 

1 + C/i 

From this last equation, taking a = 2, we get Eq. ()100|) . Substituting this result in Eq. ()101|) 
and keeping only the leading terms for d —>■ +oo, we get 

S',{u) 2C 



{d + u){l + Cfi)' 



(108) 



This equation is easily integrated and leads to the first approximation of 5*2 (x) in the large d 
limit 

S2{x) = ^—^ (109) 



{d + x2)i+Cm 

where A is an integration constant which cannot be determined explicitly at this order. We are 
now able to obtain the next order correction of a. Let us write 



a 



4C I Qi 
l + Cfi d' 



;iio) 



Inserting this expression in Eq. (jl06|) . considering the limit d — > +oo and using Eqs. (j93p and 
(jl09j) . we finally obtain 

8C(2CV-C/i-i) 



(1 + w 

This leads to the approximate expression of a to order 1/d. 

4C 



a 



1 + C/i 



2(2CV-C/i-l) 

rf(i + c^^Y 



;iii) 



112) 



This expression is valid for arbitrary values of and C, such that a <2. 
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Figure 15: The resolution of the time-dependent equations (fHH|l shows that the evolution is 
self-similar. We fix for the simulation: fi = 1, d = 3, T = 0.2, C = 0.5 and Mi = M2 = 0.5. 
For t tcoii, the rescaled densities converge to the invariant profiles Si{x) and S2{x) predicted 
by the theory. The profile of species 1 is the same as in the one-component problem : Si{x) = 



So{x) 

S2ix) 



X 



'. The profile of species 2 has been obtained by solving Eq. ()94|) numerically : 



with a = 1.66554193. 



4.3 Perturbation expansion for C ~ 1 ^^id /i ~ 1 

We now consider the case of weak asymmetry ~ 1 and C ~ 1 between the two species for any 
dimension d. In that case, 5*2 (x) will be close to Sq{x) and a will be close to 2. We set 



C = l — e, fi = 1 + 1] , a = 2 — ea^ — rja^, 



;ii3) 



S2ix) = Soix) [1 + tg^{x) + r]g^{x)] , (114) 

with €,7] <^ 1. Substituting this expansion in Eq. (jHH), it is found that the functions g({x) and 
gij,{x) satisfy the first order differential equations (for their derivatives) 

, fd + 1 \ , 2(d-2) , , 

We shall discuss these equations separately. 
4.3.1 The case /i = 1 

We first consider the case /i = 1. Equations ()85|) have been solved numerically for ^ = 1/2 
and the corresponding scaling profiles are plotted in Fig. El The numerical results lead to the 
predicted exponents : at large x, Si{x) ~ and S2{x) ~ where a is calculated using 
Eq. (jni}. We now consider the weak asymmetry limit ( = 1 — e with e ^ 1 for /i = 1 (the 
condition a < 2 imposes e > 0). Then, S2{x) = So{x)[l + e(7((x)] where 5',^(x) is the solution of 
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Eq. (jll5|) . This equation can be integrated once leading to 



The integration constant has been determined so as to satisfy the boundary condition 5'^(0) = 0. 
Now, the condition that g'(^{x) ^ as x +00, leads to an exact expression of a^. As 
^~(d+i)^x'^/2 _^ _|_^ Jqj, ^ _^ +00, the integral in Eq. pi7|) has to vanish at large x. This yields 

«c(rf) = ' L > 0. (118) 

Note that the integrals can be expressed in terms of F functions. Rewriting Eq. ()117|) in the 
form 

g'cix) = y'^'e-y'/' (^T^ - "c) dy, (119) 

we derive the large x behaviors 

(7^(x)~^ , gc{x) a^lnx. (120) 

We can also carry an expansion of a^((i) in powers of d^^ in the limit d +00. Using the 
saddle point method in Eq. ()118|) around the point y = vcT+l, we obtain 

/jx 1 3 1 15197 266999 ^fl\ ,,,,, 
= ' - 2d - - 2592M5 - 311040# ^^''^ 

We can check that the first terms of this expansion reproduce those given by Eq. ()112|1 for 
H = 1 and ( = 1 ~ e. 

4.3.2 The case C = 1 

We now consider the case C = 1- Equations ()85p have been solved numerically for /i = 2 and 
the corresponding scaling profiles are plotted in Fig. ^| They converge to the invariant profiles 
predicted by theory. We now consider the weak asymmetry limit ^ = 1 + r] with rj <^ 1 for 
C = 1 (the condition a < 2 imposes t] > 0). Then, S2{x) = So{x)[l + rig^{x)] where g^{x) 
is solution of Eq. ()116|) . Following a procedure similar to that exposed previously, we get the 
following expression of a^, 

a,id) = ^° ^ " ^ > 0. (122) 



and the asymptotic behaviors 



g'uix) , g^{x) ^ ai^lnx. (123) 



X 

The large d expansion of Eq. p22j) is 

1 5 9 23 f I \ 

and the first terms of this expansion reproduce those of Eq. ()112|) ioi ( = 1 and /i = l + rj. The 
exact values of a^((i) and a^((i) along with their 0{d'^) expansions are plotted in Fig. [T7| For 
d = 3, the exact values are a^(3) = 0.437119695 and a^(3) = 0.940162135. 
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Figure 16: The resolution of the time- dependent equations (jHH|) shows that the evolution is 
self-similar. We fix for the simulation : ( = 1, d = 3, T = 0.2, fi = 2 and Mi = M2 = 0.5. 
For t tcoii, the rescaled densities converge to the invariant profiles Si{x) and 5*2 (x) predicted 
by the theory. The profile of species 1 is the same as in the one-component problem : Si{x) = 
So{x) ~ The profile of species 2 has been obtained by solving Eq. ()94|) numerically : 

5*2 (x) ~ with a = 1.35191432. Note that the large d expansion Eq. ()l()7j) leads to a = 4/3 
in fair agreement with the exact numerical result. 




Figure 17: Numerical calculation of and given by Eqs. (|118p and (jl22j) as a function 
of the dimension d. The dashed line represents the asymptotic value for d +00. The (o) 
symbols represent the large d expansion of and to order d~'^ given in Eqs. ()121|) and 
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4.4 Other perturbation expansions 

We now consider perturbation expansions of Eq. for small and large values of /x and (. For 
/X — > and C < Cc, Eq- reduces to 

x{CSo{x) - l)52(x) + (dCSoix) - a)S2{x) = 0. (125) 

Considering the value a; = 0, we get d(So{0) — a = leading to 



Then, the scaling profile is given by 



-^^l^) = (^_2_4C + a;2)2'iC/(<i-2)- (127) 

We now wish to examine the limit [i +oo. We assume that a ^ 1/ ^ and check this scaling 
a posteriori. Using the fact that 5*2 ~ for x +oo and comparing terms of order 
in Eq. dnU, we find that 

a = -. (128) 



We note that this expression is independent on C,. The scaling profile obtained from Eq. (|94jl 
can be expressed in terms of hypergeometric functions. 
In the limit ^ — >^ 0, Eq. (jHH) simplifies into 

-xS'^i^x) + {dCS^i^x) - a)S2{x) = 0. (129) 

Considering the value x = 0, we get 

Then, the scaling profile 5*2 (x) takes the form 

= (rf_2 + a;2)2^C/(^-2)- (131) 

We now wish to examine the limit ( ^ +oo and /i > 2. Using the fact that 5*2 ~ for 
x — >• +CX) and comparing terms of order we find that /la^ — (4 + dfi)a + Ad = leading 

to a = d or a = 4/ fj,. Since a < 2, we get 

a = -. (132) 

Figure UHl shows the functions a{( = 1, /i) and a(C, = 1) obtained by solving Eq. (jMj) in d = 3 
and compares these numerical results with the asymptotic expansions obtained previously. For 
( and /i close to 1, the slope of the function a{(, fi) is given by Eqs. (|122p and (|118|) . This figure 
also confirms the asymptotic expressions ()126|) and ()13U|) obtained for /x ^ +cxd and C 0. 

5 Conclusion 

In this paper, we have extended previous studies on the thermodynamics of self-gravitating 
particles in d-dimensions to the case of multi-components systems. Our static study applies 
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Figure 18: The exponent a is plotted for /x = 1, ^ < 1 and for /i > 1, = 1 in = 3. The 
dashed hnes give the different asymptotic behaviors obtained analytically. Finally, the dotted 
line for ( = 1 corresponds to the result Eq. p07|) of the large d expansion, and is in excellent 
agreement with the exact value of a (the two curves are almost indistinguishable). 

both to the microcanonical (fixed E) and canonical (fixed T) ensembles. Thus, it describes 
ordinary stellar systems (like globular clusters) self-gravitating Brownian particles |1 and 
bacterial populations ^7||nj. We have investigated how the critical energy (Antonov point) and 
the critical temperature (Jeans point) depend on the parameters. If we take as a reference a 
single-species system with particles of mass m2 and add particles of mass mi (while removing 
some particles of mass m2 so as to keep the total mass M fixed), we find that the critical 
temperature is increased if mi > m2 and decreased if mi < m2 (an analytical estimate of the 
critical temperature has been obtained in ci > 2 using the Jeans swindle). By contrast, the 
critical energy is always increased with respect to the single species case in > 3. For given 
ratio fi = mi/m2, it presents a maximum at a certain value of x = M1/M2 ^ 0.7 (Fig. Elb). 
This maximum energy increases roughly linearly with fi (Fig. EI). As in the one-component 
case, two-dimensional systems require a specific attention. In ci = 2, there is no collapse in 
the microcanonical ensemble but there is a collapse in the canonical ensemble below a critical 
temperature. We have obtained this critical temperature analytically. For fi < 2 and for {fi > 2 
and X < X* = ~ 2)), the two species of particles form a Dirac peak at T = Tc and the 
expression (fTOj) of the critical temperature can be obtained from the Virial theorem. For (/i > 2 
and X > X* = l/(/^ ~ 2)), only the heaviest particles form a Dirac at T = Tc and the expression 
of the critical temperature (fTB)) is different. 

We have also studied the dynamics of self-gravitating Brownian particles (and bacterial pop- 
ulations) in the framework of the two-species Smoluchowski-Poisson system. This corresponds 
to the canonical ensemble. For T < Tc, there is no equilibrium state and the system collapses. 
Looking for self-similar solutions, we have shown that one species dominates the other and 
collapses as in the single-species problem with a scaling profile p(r) ~ r~^. The selection of 
the dominant species is non-trivial. For /i > 2, the dominant species is always the one with 
the heaviest particles. For /U < 2, the selection depends on the ratio ( of friction parameters as 
shown in Fig. [121 (for C = 1; the species with heaviest particles always dominates the collapse). 
The scaling profile of the "slaved" species decays with an exponent a < 2 depending on d, /i 
and X- This exponent can be calculated numerically by solving Eq. (j94j) . We have also given 
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several asymptotic expansions, see Eqs. (ITT^ . (ITTHD . (IT^ . (IT^ . (IT^ . dnUD and (ITia . 

The generalization of our approach to a continuous spectrum of masses and friction coef- 
ficients does not look straightforward. Let us focus on the simpler case of identical friction 
coefficients. The precise form of the mass spectrum is certainly highly relevant. In partic- 
ular, we expect that the value of the minimum and maximum masses (possibly and +00) 
is crucial. In addition, the behavior of the distribution near the largest mass (for instance 
p{m) ~ (m — mmax)"'*^ if ^max is finite (7 < 1), p{m) ~ otherwise (7 > 1)) is certainly an 
important ingredient. However, in the case of a bounded distribution of mass without extrav- 
agant singularities, we expect that the results obtained in the present paper will qualitatively 
hold: the heaviest particles will scale as in the one species case, while lighter particles will scale 
with a mass dependent exponent less than 2. 



A Derivation of the mean-field equations 

In this Appendix, we show that the mean-field approximation used in our study is exact in 
a proper thermodynamic limit (see |H] in the single species case). We consider the case of 
self-gravitating Brownian particles described by the stochastic equations (fTTj) . The proper 
statistical ensemble for this system is the canonical ensemble. At equilibrium, the iV-body 
distribution function is given by 

P;v(ri,vi,...,r^,v^) = ^e-'3^('-i'^^--^-^-), (133) 



where is the partition function (normalization constant) and H is the Hamiltonian 

N 

2 



1 ^ 

H = - niiv'^ + mimju{Yi — r^) = K + U, (134) 

1=1 i<j 

where u(rj — r^) = Uij = —G/[{d — 2)\ri — rj\^'^~'^^] is the gravitational potential. From Eq. ()134|) . 
it is clear that the velocity distribution is Gaussian. We shall therefore restrict ourselves to the 
configurational part 

P^(ri,...,r^) = ie-^^(^^'-''--). (135) 
We introduce the density probability for particle i of species a to be at r,, namely 

P^"\^^) = j P7v(ri,...,r;v) \{dv,. (136) 

Similarly, we define the density probability to find particle i of species a at and particle j of 
species a' at r^, 

,) = j Pjv(ri,...,riv) n ^^'^^ (137) 

The total density of particles at r is given by p(r) = ^jmj(5(r — r^). Its mean value (p(r)) = 
Y,- J mi6 (r - Ti) PAr(ri, tn) Ylj drj can be written 

X „ X X 

(pW) = EE / ^«5(r-r.)p("^(r,)dr, = 5^iV,m,p(")(r) = ^(p,(r)), (138) 

a=l i£la a=l a=l 
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where we have defined = [ga~iN + 1, QaN] with (70 = and = 1 as the interval of indices 
labehng particles of species a. In the following, we shall work with the mean density. Therefore, 
we drop the brackets (•) in order to simplify the notations. Taking the derivative of Eq. 
with respect to r,, we get 



Assuming that i & la and integrating over the other variables, we find that 



f^Y^Na' j Pt'\r,, r,)mama'^ dr,. (140) 



dr, 



Similarly, we can write an equation for Pg"" ''(rj, Yj) by integrating Eq. ()133p over A^— 2 variables. 
Then writing 

Pf "'V., r,) = P("V.)^^V.) + ^2^""'V. r,), (141) 

we can show (see jHj in the single species case) that the cumulating function P2 is of order A*""^ 
in the limit 

Aq, — > +00 with faxed Tia = — T^T^ — and Ua = . (142) 

P"~"' mx 

Thus, in this proper thermodynamic limit, we can make the mean-field approximation 

Pi""')(r„r,) = p(")(r,)p(°')(r,), (143) 

which consists in neglecting the correlations and replacing the two-body distribution function 
by a product of two one-body distribution functions. Inserting this expression in Eq. ()140|1 . we 
get 

^ = -/5m.p("V.) E / ^-'A^" \r2)ma>^dv,. (144) 

a' 

Introducing the mean density of each species 

p„(r) = Ar,m,p(")(r), (145) 

and the gravitational potential 

$(r) = E / P"(r')w(r - r) dv' , (146) 
we can rewrite the above equation in the form 

^ = -/?m,p,(r)V<l>(r). (147) 

After integration, we obtain the Boltzmann distribution 

Pa{v) = v4,e-^'"'^*('-). (148) 
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The mean potential energy W = {U) is given by 
1 ^1 

if^j ^ k=l a ^ 

imaTria' / ^(ri, r2)Mi2 (iri (ir2. (149) 

Implementing the mean-field approximation ()14Hj) . valid in the thermodynamic limit, the above 
expression simplifies into 

W = ^^iVA,m,m,, /"p(")(ri)p("')(r2)Mi2rfrit/r2, (150) 
which can be finally rewritten as 

= ly p(r)<l>(r)dr. (151) 

We now consider the dynamical problem defined by the stochastic equations ()lip. Using 
the Kramers-Moyal expansion, the Fokker-Planck equation for the evolution of the A^-body 
distribution function P/v(ri, Vi, r^r, v^r, t) reads 

^ + E (v. ■ ^ + f < . — j = — . (^A^ + 6P«v. j , (152) 

1=1 ^ ' 1=1 ^ ' 

where Fj = — (l/m.j)Vif/ is the force by unit of mass (acceleration) acting on particle i. We 
note that the stationary solution of Eq. ()152|) is the canonical distribution ()133p provided that 
the coefficients of friction and diffusion are related to each other according to the Einstein 
formula 

A = (153) 
Taking i & la and integrating over the other variables, we get 



Now 



/ = /" FiPAT Y[ drk dvk = - I {Na - l)ma^^Pt"'^ (r^, Vj, r2, V2, t) dr2 d\2 

J Na'ma'^^Pt"\ri,^i^T^2,^2,t)dr2dv2. (155) 



From the A^-body Fokker-Planck equation p52|) . we can obtain an equation for the time evo- 
lution of the two-body distribution function P2°" ""(rj, Vj, r^, Vj, t) and again show that, in the 
proper thermodynamic limit, the mean-field approximation ()143|) becomes exact. In that case, 
the expression of I simplifies into 



/ = -Pf Hr„ v„t) I 5^iV„,m„,P^"')(r2,V2,t)^rfr2rfv2 = P.^^^r., v., t)(F), 



(156) 



35 



where (F)^ = — Vj$ is the mean force (by unit of mass) acting on particle i. Introducing the 
distribution function 

/„(r, v,t) = iV„m„Pi^"^(r,v,t), (157) 
the mean-field Fokker-Planck equation takes the form 

+ v-^ + (F)-^ = Tr--U^-lS:7 + ^"/-v • (158) 



dt dr dv dv \ dv 

In the strong friction limit, the stochastic equations of motion are given by Eq. ()15|1 . In 
that case, the A^-body Fokker-Planck equation reads 

i=l ^ ' 

We note that the stationary solution of this equation is given by the configurational part of the 
canonical distribution ()133p provided that the diffusion coefficient and the mobility are related 
to each other by the Einstein relation 

D\ = ^. (160) 
Assuming that i & la and integrating over the other variables, we get 

' -^■\D'^^ + f,a I PNV.U\\dr,\. (161) 



dt dvi \ ° (9r, 



Evaluating the last term in the mean-field approximation as done previously, we find that 
which is clearly the same as 

^-V-f^Vp. + p.V$). (163) 



dt \ m, 

B Estimate of the critical temperature using the Jeans 
swindle 

In this Appendix, we extend the original Jeans instability criterion to the case of a multi- 
components system. We make the Jeans swindle, assuming that the unperturbed state is infinite 
and homogeneous. Then, we use this criterion to obtain an estimate of the critical temperature 
Tc of an inhomogeneous isothermal multi-components self-gravitating system confined within a 
box. 

Let us consider a small perturbation around an equilibrium state of the two-components 
Smoluchowski-Poisson system. The linearized equations for the perturbation can be written 

^ = iv-f^V5p, + p,V5$ + 5piV$ 

: (164) 

^^^^ ^ V ■ ( ^V5p2 + P2 V(5<l> + (5p2V$ 



dt ^2 V "^2 
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where p and $ refer to the equihbrium state. They have to be completed with the hnearized 
Poisson equation 



= SdG5p. (165) 

These equations are exact but they remain comphcated if the static solution is inhomogeneous. 
They can be solved (semi-analytically) for a one-component system [ISJ but the generalization to 
a multi-components system is not straightforward. We shall invoke here the Jeans swindle and 
consider that the equilibrium state is infinite and homogeneous although this does not rigorously 
satisfy the equations at zeroth order. With this simplifying assumption, using Eq. (jl65|) . the 
linearized equations ()164|) take the form 



85 pi 1 f ksT 

-Adpi + piSdGdp 



A5p2 + p2SdG6p . 



dt ^1 \ nil 

dSp2 1 fksT 



dt ^2 V "^2 

Writing the perturbation as 6pa ~ e*(''-'-^*), we get 



(166) 



k T \ 
-iuii + -^/c^ - SdGpi 6pi - SdGpi5p2 = 0, 
mi J 

-SdGp25pi + (-ii^^2 + -^/^^ - SdGp2] 5p2 = 0. 

V "^2 J 



(167) 



The cancellation of the determinant of the above system determines the dispersion relation. 
One can show that 7 = —iu is real so it represents either the growth rate of the perturbation 
(7 > 0) or its exponential damping (7 < 0) ^. The point of marginal stability u; = is obtained 
for the Jeans wavevector 



^2 _ SdG {niipi + m2P2) 



More generally, for a multicomponent system, we have found that 

k] = 1^ rriaPa- (169) 



ksT ^ 



The criterion of instability —iu > is equivalent to k < kj. In terms of the wavelength 
A = 27r//c, it can be written 

A > = V 170 

bdG [niipi + m2p2) 

This criterion means that if the size of the perturbation A is larger than the critical value Aj, 
the gravitational attraction will prevail over diffusion and the system will collapse. If we now 
return to our original problem, namely an isothermal system enclosed within a box of radius 
i?, a naive application of the above criterion indicates that the system is unstable if i? > Aj. 



"'^Note that if we start from the two-components barotropic Euler equations (which are the usual equations 
used in Jeans analysis) instead of the two-components Smoluchowski equations, we get the same form of dis- 
persion relation except that iLoS^a is replaced by uj^ and kBT/nia is replaced by the velocity of sound c^. In 
that case is real ; for uj'^ > 0, the system is stable and the perturbation oscillates with pulsation uj and for 
oj^ < the system is unstable and the growth rate is \iuj\. 
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Introducing the total mass of each species through the relation ~ M^/R'^, this criterion can 
be rewritten in terms of the temperature as 



As noticed in the critical temperature Tc marking the gravitational instability of box- 
confined isothermal systems can be related to the Jeans instability criterion by the above 
argument. Of course, this naive approach cannot catch the numerical constant Kd which 
appears in the expression of the critical temperature. However, this constant can be obtained 
from the numerical study of the single species case in = 3 where we have = GMm/2.52kBR 
|13j . Thus, we take = 1/2.52. Now, the interest of our treatment for a multi-components 
system is that we can obtain the dependence of the critical temperature with /x and x- Returning 
to dimensionless variables, we get the instability criterion 



where rjj is the critical inverse temperature obtained by using the Jeans swindle. This expression 
returns the single-species result for x = and for /i = 1. It is also consistent with the single- 
species result for x ~^ +oo if we redefine rj with mi instead of 1712. If we consider the limit 
/i ^ or +00 with fixed N1/N2, then x = {Ni/N2)fi — or +00, and we again recover the 
single species case. More generally, we see in Fig. 0] that this approximate expression gives a 
fair agreement with the exact solution. This is quite satisfying in view of the approximations 
made to arrive at Eq. (jl72p (we have assumed that the system is homogeneous). The relative 
success of this naive approach is explained by the fact that in d = 3 the system is weakly 
inhomogeneous at Tc. By contrast, the expression p72j) does not work at all in d = 2 (compare 
with the exact values (fTOj) and (fTBj)) because the system tends to a Dirac peak at T = Tc. 
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